目次

このレポートでは、HSP_COVID19Stress Projectの分析経過を報告します。分析の構成は以下のとおりです。分析の再現性を担保するために用いたコードも記しています。

  • (1)データハンドリングと事前分析
  • (2)階層的重回帰分析
  • (3)媒介分析

(1)データハンドリングと事前分析

1-1. ローデータの読み込み

#tidyverseパッケージ読み込み
library(tidyverse)
## -- Attaching packages ------------------------------------------ tidyverse 1.3.0 --
## √ ggplot2 3.3.0     √ purrr   0.3.3
## √ tibble  2.1.3     √ dplyr   0.8.5
## √ tidyr   1.0.2     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.5.0
## -- Conflicts --------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
#データ読み込み
lowdata <- read_csv("data.csv", na = c(".")) #セルの"."を欠損値とする
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   age = col_character(),
##   gender = col_character()
## )
## See spec(...) for full column specifications.
head(lowdata)
names(lowdata)
##  [1] "ID"          "japanese"    "age"         "gender"      "hsp1"       
##  [6] "hsp2"        "hsp3"        "hsp4"        "hsp5"        "hsp6"       
## [11] "hsp7"        "hsp8"        "hsp9"        "hsp10"       "brs1"       
## [16] "brs2"        "brs3"        "brs4"        "brs5"        "brs6"       
## [21] "covid1"      "covid2"      "covid3"      "covid4"      "covid5"     
## [26] "covid6"      "covid7"      "covid8"      "covid9"      "covid10"    
## [31] "covid11"     "covid12"     "covid13"     "covid14"     "covid15"    
## [36] "covid16"     "covid17"     "covid18"     "covid19"     "covid20"    
## [41] "covid21"     "covid22"     "covid23"     "covid24"     "covid25"    
## [46] "covid26"     "covid27"     "covid28"     "covid29"     "covid30"    
## [51] "covid31"     "covid32"     "covid33"     "covid34"     "covid35"    
## [56] "covid36"     "hsp_mean"    "hsp_total"   "brs2_t"      "brs4_t"     
## [61] "brs6_t"      "brs_mean"    "brs_total"   "covid_mean"  "covid_total"

1-2. スクリーニングと下位尺度作成

  • 日本人以外のデータを削除
  • 28歳以上のデータを削除
  • 性別無回答のデータを削除
data <- lowdata %>% filter(japanese == "0", age < 28, gender != "回答しない、その他") 
data <- data %>% filter(age != 17) #17歳のデータを削除

data <- data %>% 
  dplyr::mutate(eoe = (hsp1 + hsp2 + hsp4 + hsp6 + hsp9)/5, na.rm = TRUE) %>% #EOEの平均値
  dplyr::mutate(lst = (hsp3 + hsp7 + hsp10)/3, na.rm = TRUE) %>% #LSTの平均値
  dplyr::mutate(aes = (hsp5 + hsp8)/2, na.rm = TRUE) #AESの平均値
names(data)  
##  [1] "ID"          "japanese"    "age"         "gender"      "hsp1"       
##  [6] "hsp2"        "hsp3"        "hsp4"        "hsp5"        "hsp6"       
## [11] "hsp7"        "hsp8"        "hsp9"        "hsp10"       "brs1"       
## [16] "brs2"        "brs3"        "brs4"        "brs5"        "brs6"       
## [21] "covid1"      "covid2"      "covid3"      "covid4"      "covid5"     
## [26] "covid6"      "covid7"      "covid8"      "covid9"      "covid10"    
## [31] "covid11"     "covid12"     "covid13"     "covid14"     "covid15"    
## [36] "covid16"     "covid17"     "covid18"     "covid19"     "covid20"    
## [41] "covid21"     "covid22"     "covid23"     "covid24"     "covid25"    
## [46] "covid26"     "covid27"     "covid28"     "covid29"     "covid30"    
## [51] "covid31"     "covid32"     "covid33"     "covid34"     "covid35"    
## [56] "covid36"     "hsp_mean"    "hsp_total"   "brs2_t"      "brs4_t"     
## [61] "brs6_t"      "brs_mean"    "brs_total"   "covid_mean"  "covid_total"
## [66] "eoe"         "na.rm"       "lst"         "aes"

1-3. ヒストグラム

#ageの度数分布とヒストグラム
age_count <- dplyr::count(data, age)
knitr::kable(age_count) #テーブル化
age n
18 131
19 243
20 51
21 11
22 2
23 2
24 1
ggplot(data = data, mapping = aes(x = age, fill = factor(age))) + geom_bar() #視覚化

#genderの度数分布とヒストグラム
gender_count <- dplyr::count(data, gender)
knitr::kable(gender_count) #テーブル化
gender n
1 204
2 237
ggplot(data = data, mapping = aes(x = gender, fill = factor(gender))) + geom_bar() #視覚化

#hsp1の度数分布とヒストグラム
hsp1_count <- dplyr::count(data, hsp1)
knitr::kable(hsp1_count) #テーブル化
hsp1 n
1 15
2 54
3 72
4 37
5 159
6 69
7 35
ggplot(data = data, mapping = aes(x = hsp1, fill = factor(hsp1))) + geom_bar() #視覚化

#hsp2の度数分布とヒストグラム
hsp2_count <- dplyr::count(data, hsp2)
knitr::kable(hsp2_count) #テーブル化
hsp2 n
1 11
2 47
3 76
4 51
5 133
6 75
7 48
ggplot(data = data, mapping = aes(x = hsp2, fill = factor(hsp2))) + geom_bar() #視覚化

#hsp3の度数分布とヒストグラム
hsp3_count <- dplyr::count(data, hsp3)
knitr::kable(hsp3_count) #テーブル化
hsp3 n
1 35
2 80
3 68
4 34
5 99
6 73
7 52
ggplot(data = data, mapping = aes(x = hsp3, fill = factor(hsp3))) + geom_bar() #視覚化

#hsp4の度数分布とヒストグラム
hsp4_count <- dplyr::count(data, hsp4)
knitr::kable(hsp4_count) #テーブル化
hsp4 n
1 15
2 36
3 44
4 43
5 143
6 99
7 61
ggplot(data = data, mapping = aes(x = hsp4, fill = factor(hsp4))) + geom_bar() #視覚化

#hsp5の度数分布とヒストグラム
hsp5_count <- dplyr::count(data, hsp5)
knitr::kable(hsp5_count) #テーブル化
hsp5 n
1 23
2 65
3 73
4 70
5 92
6 72
7 46
ggplot(data = data, mapping = aes(x = hsp5, fill = factor(hsp5))) + geom_bar() #視覚化

#hsp6の度数分布とヒストグラム
hsp6_count <- dplyr::count(data, hsp6)
knitr::kable(hsp6_count) #テーブル化
hsp6 n
1 8
2 47
3 52
4 43
5 100
6 105
7 86
ggplot(data = data, mapping = aes(x = hsp6, fill = factor(hsp6))) + geom_bar() #視覚化

#hsp7の度数分布とヒストグラム
hsp7_count <- dplyr::count(data, hsp7)
knitr::kable(hsp7_count) #テーブル化
hsp7 n
1 20
2 43
3 64
4 41
5 105
6 86
7 82
ggplot(data = data, mapping = aes(x = hsp7, fill = factor(hsp7))) + geom_bar() #視覚化

#hsp8の度数分布とヒストグラム
hsp8_count <- dplyr::count(data, hsp8)
knitr::kable(hsp8_count) #テーブル化
hsp8 n
1 15
2 52
3 50
4 58
5 116
6 84
7 66
ggplot(data = data, mapping = aes(x = hsp8, fill = factor(hsp8))) + geom_bar() #視覚化

#hsp9の度数分布とヒストグラム
hsp9_count <- dplyr::count(data, hsp9)
knitr::kable(hsp9_count) #テーブル化
hsp9 n
1 16
2 30
3 47
4 59
5 125
6 86
7 78
ggplot(data = data, mapping = aes(x = hsp9, fill = factor(hsp9))) + geom_bar() #視覚化

#brs1の度数分布とヒストグラム
brs1_count <- dplyr::count(data, brs1)
knitr::kable(brs1_count) #テーブル化
brs1 n
1 50
2 138
3 85
4 140
5 28
ggplot(data = data, mapping = aes(x = brs1, fill = factor(brs1))) + geom_bar() #視覚化

#brs2の度数分布とヒストグラム
brs2_count <- dplyr::count(data, brs2)
knitr::kable(brs2_count) #テーブル化
brs2 n
1 19
2 117
3 70
4 162
5 73
ggplot(data = data, mapping = aes(x = brs2, fill = factor(brs2))) + geom_bar() #視覚化

#brs3の度数分布とヒストグラム
brs3_count <- dplyr::count(data, brs3)
knitr::kable(brs3_count) #テーブル化
brs3 n
1 48
2 150
3 77
4 135
5 31
ggplot(data = data, mapping = aes(x = brs3, fill = factor(brs3))) + geom_bar() #視覚化

#brs4の度数分布とヒストグラム
brs4_count <- dplyr::count(data, brs4)
knitr::kable(brs4_count) #テーブル化
brs4 n
1 30
2 129
3 99
4 144
5 39
ggplot(data = data, mapping = aes(x = brs4, fill = factor(brs4))) + geom_bar() #視覚化

#covid1の度数分布とヒストグラム
covid1_count <- dplyr::count(data, covid1)
knitr::kable(covid1_count) #テーブル化
covid1 n
0 79
1 154
2 158
3 32
4 18
ggplot(data = data, mapping = aes(x = covid1, fill = factor(covid1))) + geom_bar() #視覚化

#covid2の度数分布とヒストグラム
covid2_count <- dplyr::count(data, covid2)
knitr::kable(covid2_count) #テーブル化
covid2 n
0 70
1 151
2 145
3 48
4 27
ggplot(data = data, mapping = aes(x = covid2, fill = factor(covid2))) + geom_bar() #視覚化

#covid3の度数分布とヒストグラム
covid3_count <- dplyr::count(data, covid3)
knitr::kable(covid3_count) #テーブル化
covid3 n
0 86
1 196
2 122
3 28
4 9
ggplot(data = data, mapping = aes(x = covid3, fill = factor(covid3))) + geom_bar() #視覚化

#covid4の度数分布とヒストグラム
covid4_count <- dplyr::count(data, covid4)
knitr::kable(covid4_count) #テーブル化
covid4 n
0 102
1 201
2 102
3 29
4 7
ggplot(data = data, mapping = aes(x = covid4, fill = factor(covid4))) + geom_bar() #視覚化

#covid5の度数分布とヒストグラム
covid5_count <- dplyr::count(data, covid5)
knitr::kable(covid5_count) #テーブル化
covid5 n
0 84
1 216
2 91
3 37
4 13
ggplot(data = data, mapping = aes(x = covid5, fill = factor(covid5))) + geom_bar() #視覚化

#covid6の度数分布とヒストグラム
covid6_count <- dplyr::count(data, covid6)
knitr::kable(covid6_count) #テーブル化
covid6 n
0 68
1 195
2 105
3 60
4 13
ggplot(data = data, mapping = aes(x = covid6, fill = factor(covid6))) + geom_bar() #視覚化

#covid7の度数分布とヒストグラム
covid7_count <- dplyr::count(data, covid7)
knitr::kable(covid7_count) #テーブル化
covid7 n
0 211
1 179
2 37
3 8
4 6
ggplot(data = data, mapping = aes(x = covid7, fill = factor(covid7))) + geom_bar() #視覚化

#covid8の度数分布とヒストグラム
covid8_count <- dplyr::count(data, covid8)
knitr::kable(covid8_count) #テーブル化
covid8 n
0 221
1 157
2 45
3 14
4 4
ggplot(data = data, mapping = aes(x = covid8, fill = factor(covid8))) + geom_bar() #視覚化

#covid9の度数分布とヒストグラム
covid9_count <- dplyr::count(data, covid9)
knitr::kable(covid9_count) #テーブル化
covid9 n
0 173
1 151
2 74
3 36
4 7
ggplot(data = data, mapping = aes(x = covid9, fill = factor(covid9))) + geom_bar() #視覚化

#covid10の度数分布とヒストグラム
covid10_count <- dplyr::count(data, covid10)
knitr::kable(covid10_count) #テーブル化
covid10 n
0 194
1 145
2 68
3 30
4 4
ggplot(data = data, mapping = aes(x = covid10, fill = factor(covid10))) + geom_bar() #視覚化

#covid11の度数分布とヒストグラム
covid11_count <- dplyr::count(data, covid11)
knitr::kable(covid11_count) #テーブル化
covid11 n
0 217
1 147
2 45
3 28
4 4
ggplot(data = data, mapping = aes(x = covid11, fill = factor(covid11))) + geom_bar() #視覚化

#covid12の度数分布とヒストグラム
covid12_count <- dplyr::count(data, covid12)
knitr::kable(covid12_count) #テーブル化
covid12 n
0 223
1 142
2 40
3 28
4 8
ggplot(data = data, mapping = aes(x = covid12, fill = factor(covid12))) + geom_bar() #視覚化

#covid13の度数分布とヒストグラム
covid13_count <- dplyr::count(data, covid13)
knitr::kable(covid13_count) #テーブル化
covid13 n
0 168
1 152
2 71
3 38
4 12
ggplot(data = data, mapping = aes(x = covid13, fill = factor(covid13))) + geom_bar() #視覚化

#covid14の度数分布とヒストグラム
covid14_count <- dplyr::count(data, covid14)
knitr::kable(covid14_count) #テーブル化
covid14 n
0 241
1 135
2 40
3 22
4 3
ggplot(data = data, mapping = aes(x = covid14, fill = factor(covid14))) + geom_bar() #視覚化

#covid15の度数分布とヒストグラム
covid15_count <- dplyr::count(data, covid15)
knitr::kable(covid15_count) #テーブル化
covid15 n
0 223
1 122
2 64
3 28
4 4
ggplot(data = data, mapping = aes(x = covid15, fill = factor(covid15))) + geom_bar() #視覚化

#covid16の度数分布とヒストグラム
covid16_count <- dplyr::count(data, covid16)
knitr::kable(covid16_count) #テーブル化
covid16 n
0 153
1 142
2 71
3 55
4 20
ggplot(data = data, mapping = aes(x = covid16, fill = factor(covid16))) + geom_bar() #視覚化

#covid17の度数分布とヒストグラム
covid17_count <- dplyr::count(data, covid17)
knitr::kable(covid17_count) #テーブル化
covid17 n
0 161
1 138
2 76
3 50
4 16
ggplot(data = data, mapping = aes(x = covid17, fill = factor(covid17))) + geom_bar() #視覚化

#covid18の度数分布とヒストグラム
covid18_count <- dplyr::count(data, covid18)
knitr::kable(covid18_count) #テーブル化
covid18 n
0 223
1 138
2 52
3 21
4 7
ggplot(data = data, mapping = aes(x = covid18, fill = factor(covid18))) + geom_bar() #視覚化

#covid19の度数分布とヒストグラム
covid19_count <- dplyr::count(data, covid19)
knitr::kable(covid19_count) #テーブル化
covid19 n
0 48
1 107
2 140
3 94
4 52
ggplot(data = data, mapping = aes(x = covid19, fill = factor(covid19))) + geom_bar() #視覚化

#covid20の度数分布とヒストグラム
covid20_count <- dplyr::count(data, covid20)
knitr::kable(covid20_count) #テーブル化
covid20 n
0 47
1 98
2 134
3 113
4 49
ggplot(data = data, mapping = aes(x = covid20, fill = factor(covid20))) + geom_bar() #視覚化

#covid21の度数分布とヒストグラム
covid21_count <- dplyr::count(data, covid21)
knitr::kable(covid21_count) #テーブル化
covid21 n
0 57
1 115
2 142
3 95
4 32
ggplot(data = data, mapping = aes(x = covid21, fill = factor(covid21))) + geom_bar() #視覚化

#covid22の度数分布とヒストグラム
covid22_count <- dplyr::count(data, covid22)
knitr::kable(covid22_count) #テーブル化
covid22 n
0 159
1 177
2 60
3 28
4 17
ggplot(data = data, mapping = aes(x = covid22, fill = factor(covid22))) + geom_bar() #視覚化

#covid23の度数分布とヒストグラム
covid23_count <- dplyr::count(data, covid23)
knitr::kable(covid23_count) #テーブル化
covid23 n
0 163
1 174
2 60
3 34
4 10
ggplot(data = data, mapping = aes(x = covid23, fill = factor(covid23))) + geom_bar() #視覚化

#covid24の度数分布とヒストグラム
covid24_count <- dplyr::count(data, covid24)
knitr::kable(covid24_count) #テーブル化
covid24 n
0 174
1 163
2 65
3 30
4 9
ggplot(data = data, mapping = aes(x = covid24, fill = factor(covid24))) + geom_bar() #視覚化

#covid25の度数分布とヒストグラム
covid25_count <- dplyr::count(data, covid25)
knitr::kable(covid25_count) #テーブル化
covid25 n
0 274
1 132
2 24
3 8
4 3
ggplot(data = data, mapping = aes(x = covid25, fill = factor(covid25))) + geom_bar() #視覚化

#covid26の度数分布とヒストグラム
covid26_count <- dplyr::count(data, covid26)
knitr::kable(covid26_count) #テーブル化
covid26 n
0 212
1 135
2 57
3 28
4 9
ggplot(data = data, mapping = aes(x = covid26, fill = factor(covid26))) + geom_bar() #視覚化

#covid27の度数分布とヒストグラム
covid27_count <- dplyr::count(data, covid27)
knitr::kable(covid27_count) #テーブル化
covid27 n
0 365
1 64
2 8
3 3
4 1
ggplot(data = data, mapping = aes(x = covid27, fill = factor(covid27))) + geom_bar() #視覚化

#covid28の度数分布とヒストグラム
covid28_count <- dplyr::count(data, covid28)
knitr::kable(covid28_count) #テーブル化
covid28 n
0 274
1 95
2 45
3 18
4 9
ggplot(data = data, mapping = aes(x = covid28, fill = factor(covid28))) + geom_bar() #視覚化

#covid29の度数分布とヒストグラム
covid29_count <- dplyr::count(data, covid29)
knitr::kable(covid29_count) #テーブル化
covid29 n
0 352
1 65
2 15
3 7
4 2
ggplot(data = data, mapping = aes(x = covid29, fill = factor(covid29))) + geom_bar() #視覚化

#covid30の度数分布とヒストグラム
covid30_count <- dplyr::count(data, covid30)
knitr::kable(covid30_count) #テーブル化
covid30 n
0 372
1 49
2 10
3 6
4 4
ggplot(data = data, mapping = aes(x = covid30, fill = factor(covid30))) + geom_bar() #視覚化

#covid31の度数分布とヒストグラム
covid31_count <- dplyr::count(data, covid31)
knitr::kable(covid31_count) #テーブル化
covid31 n
0 230
1 110
2 70
3 25
4 6
ggplot(data = data, mapping = aes(x = covid31, fill = factor(covid31))) + geom_bar() #視覚化

#covid32の度数分布とヒストグラム
covid32_count <- dplyr::count(data, covid32)
knitr::kable(covid32_count) #テーブル化
covid32 n
0 357
1 69
2 11
3 2
4 2
ggplot(data = data, mapping = aes(x = covid32, fill = factor(covid32))) + geom_bar() #視覚化

#covid33の度数分布とヒストグラム
covid33_count <- dplyr::count(data, covid33)
knitr::kable(covid33_count) #テーブル化
covid33 n
0 277
1 83
2 51
3 22
4 8
ggplot(data = data, mapping = aes(x = covid33, fill = factor(covid33))) + geom_bar() #視覚化

#covid34の度数分布とヒストグラム
covid34_count <- dplyr::count(data, covid34)
knitr::kable(covid34_count) #テーブル化
covid34 n
0 134
1 69
2 104
3 85
4 49
ggplot(data = data, mapping = aes(x = covid34, fill = factor(covid34))) + geom_bar() #視覚化

#covid35の度数分布とヒストグラム
covid35_count <- dplyr::count(data, covid35)
knitr::kable(covid35_count) #テーブル化
covid35 n
0 232
1 96
2 61
3 35
4 17
ggplot(data = data, mapping = aes(x = covid35, fill = factor(covid35))) + geom_bar() #視覚化

#covid36の度数分布とヒストグラム
covid36_count <- dplyr::count(data, covid36)
knitr::kable(covid36_count) #テーブル化
covid36 n
0 172
1 75
2 88
3 69
4 37
ggplot(data = data, mapping = aes(x = covid36, fill = factor(covid36))) + geom_bar() #視覚化

#hsp_meanの度数分布とヒストグラム
hsp_mean_count <- dplyr::count(data, hsp_mean)
knitr::kable(hsp_mean_count) #テーブル化
hsp_mean n
1.2 1
1.3 1
1.7 1
1.8 1
2.0 2
2.3 6
2.4 1
2.5 4
2.6 2
2.7 8
2.8 7
2.9 8
3.0 4
3.1 8
3.2 5
3.3 13
3.4 10
3.5 15
3.6 9
3.7 11
3.8 16
3.9 10
4.0 12
4.1 12
4.2 16
4.3 15
4.4 10
4.5 18
4.6 16
4.7 10
4.8 19
4.9 14
5.0 17
5.1 17
5.2 20
5.3 9
5.4 14
5.5 4
5.6 16
5.7 6
5.8 6
5.9 3
6.0 5
6.1 5
6.2 6
6.3 4
6.4 2
6.5 3
6.6 6
6.7 1
6.8 11
7.0 1
ggplot(data = data, mapping = aes(x = hsp_mean, fill = factor(hsp_mean))) + geom_histogram(binwidth = 0.5) + guides(fill = "none") #視覚化

#eoeの度数分布とヒストグラム
eoe_mean_count <- dplyr::count(data, eoe)
knitr::kable(eoe_mean_count) #テーブル化
eoe n
1.0 1
1.4 2
1.8 1
2.0 5
2.2 3
2.4 8
2.6 9
2.8 5
3.0 14
3.2 12
3.4 14
3.6 16
3.8 19
4.0 23
4.2 22
4.4 32
4.6 22
4.8 23
5.0 33
5.2 37
5.4 26
5.6 13
5.8 24
6.0 22
6.2 17
6.4 16
6.6 12
6.8 3
7.0 7
ggplot(data = data, mapping = aes(x = eoe, fill = factor(eoe))) + geom_histogram(binwidth = 0.5) + guides(fill = "none") #視覚化

#lstの度数分布とヒストグラム
lst_mean_count <- dplyr::count(data, lst)
knitr::kable(lst_mean_count) #テーブル化
lst n
1.000000 12
1.333333 5
1.666667 15
2.000000 21
2.333333 29
2.666667 28
3.000000 28
3.333333 19
3.666667 27
4.000000 34
4.333333 30
4.666667 41
5.000000 23
5.333333 32
5.666667 21
6.000000 28
6.333333 12
6.666667 10
7.000000 26
ggplot(data = data, mapping = aes(x = lst, fill = factor(lst))) + geom_histogram(binwidth = 0.5) + guides(fill = "none") #視覚化

#aesの度数分布とヒストグラム
aes_mean_count <- dplyr::count(data, aes)
knitr::kable(aes_mean_count) #テーブル化
aes n
1.0 9
1.5 8
2.0 28
2.5 25
3.0 33
3.5 44
4.0 48
4.5 54
5.0 47
5.5 51
6.0 32
6.5 25
7.0 37
ggplot(data = data, mapping = aes(x = aes, fill = factor(aes))) + geom_histogram(binwidth = 0.5) + guides(fill = "none") #視覚化

#brs_meanの度数分布とヒストグラム
brs_mean_count <- dplyr::count(data, brs_mean)
knitr::kable(brs_mean_count) #テーブル化
brs_mean n
1.000000 6
1.166667 4
1.333333 6
1.500000 12
1.666667 7
1.833333 15
2.000000 22
2.166667 32
2.333333 23
2.500000 35
2.666667 28
2.833333 25
3.000000 33
3.166667 30
3.333333 15
3.500000 30
3.666667 27
3.833333 22
4.000000 30
4.166667 19
4.333333 9
4.500000 5
4.666667 3
4.833333 1
5.000000 2
ggplot(data = data, mapping = aes(x = brs_mean, fill = factor(brs_mean))) + geom_histogram(binwidth = 0.5) + guides(fill = "none") #視覚化

#covid_meanの度数分布とヒストグラム
covid_mean_count <- dplyr::count(data, covid_mean)
knitr::kable(covid_mean_count) #テーブル化
covid_mean n
0.0000000 3
0.0277778 5
0.0555556 2
0.0833333 1
0.1111111 2
0.1388889 1
0.1666667 1
0.1944444 4
0.2222222 3
0.2500000 5
0.2777778 4
0.3055556 2
0.3333333 7
0.3611111 7
0.3888889 4
0.4166667 5
0.4444444 12
0.4722222 6
0.5000000 4
0.5277778 8
0.5555556 7
0.5833333 7
0.6111111 13
0.6388889 8
0.6666667 9
0.6944444 16
0.7222222 10
0.7500000 7
0.7777778 10
0.8055556 10
0.8333333 8
0.8611111 6
0.8888889 10
0.9166667 10
0.9444444 10
0.9722222 12
1.0000000 12
1.0277778 5
1.0555556 10
1.0833333 9
1.1111111 6
1.1388889 4
1.1666667 12
1.1944444 5
1.2222222 8
1.2500000 7
1.2777778 5
1.3055556 10
1.3333333 7
1.3611111 5
1.3888889 8
1.4166667 7
1.4444444 8
1.4722222 6
1.5000000 4
1.5277778 4
1.5833333 2
1.6111111 1
1.6388889 3
1.6666667 6
1.6944444 6
1.7222222 4
1.7500000 2
1.7777778 3
1.8055556 1
1.8333333 1
1.8611111 3
1.8888889 4
1.9166667 3
1.9444444 3
1.9722222 4
2.0277778 3
2.0833333 1
2.1388889 1
2.2222222 3
2.2777778 1
2.3888889 1
2.4166667 1
2.6388889 1
3.0000000 1
3.0555556 1
ggplot(data = data, mapping = aes(x = covid_mean, fill = factor(covid_mean))) + geom_histogram(binwidth = 0.5) + guides(fill = "none") #視覚化

1-4. 記述統計量

#年齢
data$age <- as.integer(data$age)
mean(data$age) #平均年齢
## [1] 18.91156
sd(data$age) #SD
## [1] 0.8223456
min(data$age)
## [1] 18
max(data$age)
## [1] 24
#HSP、レジリエンス、コロナストレス
stat <- data %>% 
  dplyr::summarise(hsp_m = mean(hsp_mean, na.rm = TRUE),
                   hsp_sd = sd(hsp_mean, na.rm = TRUE),
                   eoe_m = mean(eoe, na.rm = TRUE),
                   eoe_sd = sd(eoe, na.rm = TRUE),
                   lst_m = mean(lst, na.rm = TRUE),
                   lst_sd = sd(lst, na.rm = TRUE),
                   aes_m = mean(aes, na.rm = TRUE),
                   aes_sd = sd(aes, na.rm = TRUE),
                   brs_m = mean(brs_mean, na.rm = TRUE),
                   brs_sd = sd(brs_mean, na.rm = TRUE),
                   covid_m = mean(covid_mean, na.rm = TRUE),
                   covid_sd = sd(covid_mean, na.rm = TRUE))
knitr::kable(stat, digits = 2) #出力
hsp_m hsp_sd eoe_m eoe_sd lst_m lst_sd aes_m aes_sd brs_m brs_sd covid_m covid_sd
4.48 1.11 4.7 1.19 4.15 1.59 4.44 1.55 2.94 0.85 0.99 0.52
min(data$hsp_mean)
## [1] 1.2
max(data$hsp_mean)
## [1] 7
min(data$brs_mean)
## [1] 1
max(data$brs_mean)
## [1] 5
min(data$covid_mean)
## [1] 0
max(data$covid_mean)
## [1] 3.055556

1-5. 内的一貫性の算出

library(psych)
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(GPArotation)

#HSP
omega(data[, c(5:14)],1,fm="ml") #alpha 0.85 #HSP
## Omega_h for 1 factor is not meaningful, just omega_t
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.85 
## G.6:                   0.87 
## Omega Hierarchical:    0.84 
## Omega H asymptotic:    0.98 
## Omega Total            0.85 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##          g  F1*   h2   u2 p2
## hsp1  0.64      0.41 0.59  1
## hsp2  0.73      0.54 0.46  1
## hsp3  0.76      0.58 0.42  1
## hsp4  0.60      0.36 0.64  1
## hsp5  0.41      0.17 0.83  1
## hsp6  0.54      0.29 0.71  1
## hsp7  0.74      0.55 0.45  1
## hsp8  0.40      0.16 0.84  1
## hsp9  0.47      0.22 0.78  1
## hsp10 0.67      0.45 0.55  1
## 
## With eigenvalues of:
##   g F1* 
## 3.7 0.0 
## 
## general/max  2.694489e+16   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 35  and the fit is  1.09 
## The number of observations was  441  with Chi Square =  476.49  with prob <  7.1e-79
## The root mean square of the residuals is  0.1 
## The df corrected root mean square of the residuals is  0.12
## RMSEA index =  0.169  and the 10 % confidence intervals are  0.156 0.183
## BIC =  263.38
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 35  and the fit is  1.09 
## The number of observations was  441  with Chi Square =  476.49  with prob <  7.1e-79
## The root mean square of the residuals is  0.1 
## The df corrected root mean square of the residuals is  0.12 
## 
## RMSEA index =  0.169  and the 10 % confidence intervals are  0.156 0.183
## BIC =  263.38 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.94   0
## Multiple R square of scores with factors      0.87   0
## Minimum correlation of factor score estimates 0.75  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.85 0.84
## Omega general for total scores and subscales  0.84 0.84
## Omega group for total scores and subscales    0.00 0.00
omega(data[, c(5,6,8,10,13)],1,fm="ml") #alpha 0.79 #EOE
## Omega_h for 1 factor is not meaningful, just omega_t
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.79 
## G.6:                   0.77 
## Omega Hierarchical:    0.79 
## Omega H asymptotic:    0.99 
## Omega Total            0.8 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##         g  F1*   h2   u2 p2
## hsp1 0.68      0.46 0.54  1
## hsp2 0.85      0.72 0.28  1
## hsp4 0.67      0.45 0.55  1
## hsp6 0.54      0.29 0.71  1
## hsp9 0.54      0.29 0.71  1
## 
## With eigenvalues of:
##   g F1* 
## 2.2 0.0 
## 
## general/max  4.006494e+16   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 5  and the fit is  0.05 
## The number of observations was  441  with Chi Square =  21.07  with prob <  0.00079
## The root mean square of the residuals is  0.04 
## The df corrected root mean square of the residuals is  0.06
## RMSEA index =  0.085  and the 10 % confidence intervals are  0.05 0.125
## BIC =  -9.37
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 5  and the fit is  0.05 
## The number of observations was  441  with Chi Square =  21.07  with prob <  0.00079
## The root mean square of the residuals is  0.04 
## The df corrected root mean square of the residuals is  0.06 
## 
## RMSEA index =  0.085  and the 10 % confidence intervals are  0.05 0.125
## BIC =  -9.37 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.91   0
## Multiple R square of scores with factors      0.84   0
## Minimum correlation of factor score estimates 0.67  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.80 0.79
## Omega general for total scores and subscales  0.79 0.79
## Omega group for total scores and subscales    0.00 0.00
omega(data[, c(7,11,14)],1,fm="ml") #alpha 0.84 #LST
## Omega_h for 1 factor is not meaningful, just omega_t
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.84 
## G.6:                   0.78 
## Omega Hierarchical:    0.84 
## Omega H asymptotic:    1 
## Omega Total            0.84 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##          g  F1*   h2   u2 p2
## hsp3  0.84      0.71 0.29  1
## hsp7  0.86      0.75 0.25  1
## hsp10 0.69      0.48 0.52  1
## 
## With eigenvalues of:
##   g F1* 
## 1.9 0.0 
## 
## general/max  3.476011e+16   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 0  and the fit is  0 
## The number of observations was  441  with Chi Square =  0  with prob <  NA
## The root mean square of the residuals is  0 
## The df corrected root mean square of the residuals is  NA
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 0  and the fit is  0 
## The number of observations was  441  with Chi Square =  0  with prob <  NA
## The root mean square of the residuals is  0 
## The df corrected root mean square of the residuals is  NA 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.93   0
## Multiple R square of scores with factors      0.86   0
## Minimum correlation of factor score estimates 0.72  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.84 0.84
## Omega general for total scores and subscales  0.84 0.84
## Omega group for total scores and subscales    0.00 0.00
omega(data[, c(9,12)],1,fm="ml") #alpha 0.78 #AES
## Omega_h for 1 factor is not meaningful, just omega_t
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.78 
## G.6:                   0.64 
## Omega Hierarchical:    0.78 
## Omega H asymptotic:    1 
## Omega Total            0.78 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##        g  F1*   h2   u2 p2
## hsp5 0.8      0.64 0.36  1
## hsp8 0.8      0.64 0.36  1
## 
## With eigenvalues of:
##   g F1* 
## 1.3 0.0 
## 
## general/max  Inf   max/min =   NaN
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are -1  and the fit is  0 
## The number of observations was  441  with Chi Square =  0  with prob <  NA
## The root mean square of the residuals is  0 
## The df corrected root mean square of the residuals is  NA
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are -1  and the fit is  0 
## The number of observations was  441  with Chi Square =  0  with prob <  NA
## The root mean square of the residuals is  0 
## The df corrected root mean square of the residuals is  NA 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.88   0
## Multiple R square of scores with factors      0.78   0
## Minimum correlation of factor score estimates 0.56  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.78 0.78
## Omega general for total scores and subscales  0.78 0.78
## Omega group for total scores and subscales    0.00 0.00
#BRS
omega(data[, c(15,59,17,60,19,61)],1,fm="ml") #alpha 0.85
## Omega_h for 1 factor is not meaningful, just omega_t
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.85 
## G.6:                   0.84 
## Omega Hierarchical:    0.85 
## Omega H asymptotic:    1 
## Omega Total            0.85 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##           g  F1*   h2   u2 p2
## brs1   0.80      0.63 0.37  1
## brs2_t 0.74      0.55 0.45  1
## brs3   0.85      0.72 0.28  1
## brs4_t 0.64      0.41 0.59  1
## brs5   0.70      0.49 0.51  1
## brs6_t 0.43      0.18 0.82  1
## 
## With eigenvalues of:
##   g F1* 
##   3   0 
## 
## general/max  2.680442e+16   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 9  and the fit is  0.11 
## The number of observations was  441  with Chi Square =  48.92  with prob <  1.7e-07
## The root mean square of the residuals is  0.04 
## The df corrected root mean square of the residuals is  0.06
## RMSEA index =  0.1  and the 10 % confidence intervals are  0.074 0.129
## BIC =  -5.88
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9  and the fit is  0.11 
## The number of observations was  441  with Chi Square =  48.92  with prob <  1.7e-07
## The root mean square of the residuals is  0.04 
## The df corrected root mean square of the residuals is  0.06 
## 
## RMSEA index =  0.1  and the 10 % confidence intervals are  0.074 0.129
## BIC =  -5.88 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.94   0
## Multiple R square of scores with factors      0.88   0
## Minimum correlation of factor score estimates 0.76  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.85 0.85
## Omega general for total scores and subscales  0.85 0.85
## Omega group for total scores and subscales    0.00 0.00
#COVID19
omega(data[, c(21:56)],1,fm="ml") #alpha 0.93
## Omega_h for 1 factor is not meaningful, just omega_t
## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.93 
## G.6:                   0.96 
## Omega Hierarchical:    0.91 
## Omega H asymptotic:    0.98 
## Omega Total            0.93 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##            g  F1*   h2   u2 p2
## covid1  0.50      0.25 0.75  1
## covid2  0.50      0.25 0.75  1
## covid3  0.48      0.23 0.77  1
## covid4  0.54      0.29 0.71  1
## covid5  0.48      0.23 0.77  1
## covid6  0.46      0.21 0.79  1
## covid7  0.43      0.18 0.82  1
## covid8  0.37      0.13 0.87  1
## covid9  0.54      0.29 0.71  1
## covid10 0.43      0.19 0.81  1
## covid11 0.47      0.23 0.77  1
## covid12 0.48      0.23 0.77  1
## covid13 0.57      0.32 0.68  1
## covid14 0.66      0.43 0.57  1
## covid15 0.69      0.48 0.52  1
## covid16 0.63      0.40 0.60  1
## covid17 0.68      0.46 0.54  1
## covid18 0.59      0.35 0.65  1
## covid19 0.56      0.31 0.69  1
## covid20 0.56      0.31 0.69  1
## covid21 0.62      0.39 0.61  1
## covid22 0.57      0.33 0.67  1
## covid23 0.60      0.35 0.65  1
## covid24 0.59      0.34 0.66  1
## covid25 0.61      0.37 0.63  1
## covid26 0.62      0.39 0.61  1
## covid27 0.48      0.23 0.77  1
## covid28 0.47      0.22 0.78  1
## covid29 0.49      0.24 0.76  1
## covid30 0.46      0.21 0.79  1
## covid31 0.38      0.14 0.86  1
## covid32 0.31      0.09 0.91  1
## covid33 0.25      0.06 0.94  1
## covid34 0.40      0.16 0.84  1
## covid35 0.51      0.26 0.74  1
## covid36 0.27      0.07 0.93  1
## 
## With eigenvalues of:
##   g F1* 
## 9.6 0.0 
## 
## general/max  1.388817e+16   max/min =   1
## mean percent general =  1    with sd =  0 and cv of  0 
## Explained Common Variance of the general factor =  1 
## 
## The degrees of freedom are 594  and the fit is  13.48 
## The number of observations was  441  with Chi Square =  5748.37  with prob <  0
## The root mean square of the residuals is  0.13 
## The df corrected root mean square of the residuals is  0.14
## RMSEA index =  0.14  and the 10 % confidence intervals are  0.137 0.144
## BIC =  2131.47
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 594  and the fit is  13.48 
## The number of observations was  441  with Chi Square =  5748.37  with prob <  0
## The root mean square of the residuals is  0.13 
## The df corrected root mean square of the residuals is  0.14 
## 
## RMSEA index =  0.14  and the 10 % confidence intervals are  0.137 0.144
## BIC =  2131.47 
## 
## Measures of factor score adequacy             
##                                                  g F1*
## Correlation of scores with factors            0.97   0
## Multiple R square of scores with factors      0.93   0
## Minimum correlation of factor score estimates 0.87  -1
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*
## Omega total for total scores and subscales    0.93 0.91
## Omega general for total scores and subscales  0.91 0.91
## Omega group for total scores and subscales    0.00 0.00

1-6. 相関係数の算出

cor.test(data$hsp_mean, data$brs_mean)#HSP-BRS
## 
##  Pearson's product-moment correlation
## 
## data:  data$hsp_mean and data$brs_mean
## t = -11.311, df = 439, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5442654 -0.3993666
## sample estimates:
##        cor 
## -0.4750297
cor.test(data$hsp_mean, data$covid_mean)#HSP-COVID
## 
##  Pearson's product-moment correlation
## 
## data:  data$hsp_mean and data$covid_mean
## t = 5.4046, df = 439, p-value = 1.068e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1601282 0.3353283
## sample estimates:
##       cor 
## 0.2497714
cor.test(data$brs_mean, data$covid_mean)#BRS-COVID
## 
##  Pearson's product-moment correlation
## 
## data:  data$brs_mean and data$covid_mean
## t = -4.86, df = 439, p-value = 1.637e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3127336 -0.1354337
## sample estimates:
##        cor 
## -0.2259541
#HSP-BRSプロット
ggplot(data, aes(x = hsp_mean, y = brs_mean, color = gender)) + geom_point() + geom_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'

#HSP-COVIDプロット
ggplot(data, aes(x = hsp_mean, y = covid_mean, color = gender)) + geom_point() + geom_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'

#BRS-COVIDプロット
ggplot(data, aes(x = brs_mean, y = covid_mean, color = gender)) + geom_point() + geom_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'

cor.test(data$eoe, data$brs_mean)#EOE-BRS
## 
##  Pearson's product-moment correlation
## 
## data:  data$eoe and data$brs_mean
## t = -14.457, df = 439, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6279877 -0.5011084
## sample estimates:
##        cor 
## -0.5679123
cor.test(data$eoe, data$covid_mean)#EOE-COVID
## 
##  Pearson's product-moment correlation
## 
## data:  data$eoe and data$covid_mean
## t = 5.1951, df = 439, p-value = 3.139e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1506684 0.3266966
## sample estimates:
##       cor 
## 0.2406604
cor.test(data$lst, data$brs_mean)#LST-BRS
## 
##  Pearson's product-moment correlation
## 
## data:  data$lst and data$brs_mean
## t = -7.346, df = 439, p-value = 1.003e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4115237 -0.2450529
## sample estimates:
##        cor 
## -0.3308598
cor.test(data$lst, data$covid_mean)#LST-COVID
## 
##  Pearson's product-moment correlation
## 
## data:  data$lst and data$covid_mean
## t = 4.6977, df = 439, p-value = 3.526e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1280126 0.3059039
## sample estimates:
##       cor 
## 0.2187753
cor.test(data$aes, data$brs_mean)#AES-BRS
## 
##  Pearson's product-moment correlation
## 
## data:  data$aes and data$brs_mean
## t = -2.1412, df = 439, p-value = 0.03281
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.193208196 -0.008366039
## sample estimates:
##        cor 
## -0.1016645
cor.test(data$aes, data$covid_mean)#AES-COVID
## 
##  Pearson's product-moment correlation
## 
## data:  data$aes and data$covid_mean
## t = 2.0213, df = 439, p-value = 0.04386
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.002669909 0.187718535
## sample estimates:
##        cor 
## 0.09602384
cor.test(data$eoe, data$lst)#EOE-LST
## 
##  Pearson's product-moment correlation
## 
## data:  data$eoe and data$lst
## t = 14.632, df = 439, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5062368 0.6321315
## sample estimates:
##       cor 
## 0.5725495
cor.test(data$eoe, data$aes)#EOE-AES
## 
##  Pearson's product-moment correlation
## 
## data:  data$eoe and data$aes
## t = 7.7999, df = 439, p-value = 4.558e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2641046 0.4283034
## sample estimates:
##       cor 
## 0.3488786
cor.test(data$lst, data$aes)#LST-AES
## 
##  Pearson's product-moment correlation
## 
## data:  data$lst and data$aes
## t = 7.6918, df = 439, p-value = 9.636e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2595966 0.4243432
## sample estimates:
##       cor 
## 0.3446207
cor.test(data$age, data$hsp_mean)#AGE-HSP
## 
##  Pearson's product-moment correlation
## 
## data:  data$age and data$hsp_mean
## t = -1.5563, df = 439, p-value = 0.1203
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.16630361  0.01943625
## sample estimates:
##         cor 
## -0.07407607
cor.test(data$age, data$brs_mean)#AGE-BRS
## 
##  Pearson's product-moment correlation
## 
## data:  data$age and data$brs_mean
## t = 1.6534, df = 439, p-value = 0.09896
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.01481838  0.17079160
## sample estimates:
##        cor 
## 0.07866834
cor.test(data$age, data$covid_mean)#AGE-COVID
## 
##  Pearson's product-moment correlation
## 
## data:  data$age and data$covid_mean
## t = -1.0962, df = 439, p-value = 0.2736
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.14491823  0.04133214
## sample estimates:
##         cor 
## -0.05224738
data$gender <- as.integer(data$gender)
cor.test(data$gender, data$hsp_mean)#GENDER-HSP
## 
##  Pearson's product-moment correlation
## 
## data:  data$gender and data$hsp_mean
## t = 3.7229, df = 439, p-value = 0.0002226
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.08292007 0.26400845
## sample estimates:
##       cor 
## 0.1749434
cor.test(data$gender, data$brs_mean)#GENDER-BRS
## 
##  Pearson's product-moment correlation
## 
## data:  data$gender and data$brs_mean
## t = -1.5585, df = 439, p-value = 0.1198
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.16640535  0.01933165
## sample estimates:
##         cor 
## -0.07418013
cor.test(data$gender, data$covid_mean)#GENDER-COVID
## 
##  Pearson's product-moment correlation
## 
## data:  data$gender and data$covid_mean
## t = 2.6364, df = 439, p-value = 0.008676
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.03183731 0.21570726
## sample estimates:
##      cor 
## 0.124844

1-7. 性差

#差の検定
t.test(hsp_mean ~ gender, data = data, var.equal = TRUE) #HSP
## 
##  Two Sample t-test
## 
## data:  hsp_mean by gender
## t = -3.7229, df = 439, p-value = 0.0002226
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.5947566 -0.1837617
## sample estimates:
## mean in group 1 mean in group 2 
##        4.271078        4.660338
t.test(eoe ~ gender, data = data, var.equal = TRUE) #EOE
## 
##  Two Sample t-test
## 
## data:  eoe by gender
## t = -3.8322, df = 439, p-value = 0.0001455
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6498287 -0.2092430
## sample estimates:
## mean in group 1 mean in group 2 
##        4.466667        4.896203
t.test(lst ~ gender, data = data, var.equal = TRUE) #LST
## 
##  Two Sample t-test
## 
## data:  lst by gender
## t = -2.4208, df = 439, p-value = 0.01589
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.66219593 -0.06875467
## sample estimates:
## mean in group 1 mean in group 2 
##        3.950980        4.316456
t.test(aes ~ gender, data = data, var.equal = TRUE) #AES
## 
##  Two Sample t-test
## 
## data:  aes by gender
## t = -2.2055, df = 439, p-value = 0.02793
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.61317953 -0.03530644
## sample estimates:
## mean in group 1 mean in group 2 
##        4.262255        4.586498
t.test(brs_mean ~ gender, data = data, var.equal = TRUE) #BRS
## 
##  Two Sample t-test
## 
## data:  brs_mean by gender
## t = 1.5585, df = 439, p-value = 0.1198
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0330620  0.2863713
## sample estimates:
## mean in group 1 mean in group 2 
##        3.010621        2.883966
t.test(covid_mean ~ gender, data = data, var.equal = TRUE) #COVID
## 
##  Two Sample t-test
## 
## data:  covid_mean by gender
## t = -2.6364, df = 439, p-value = 0.008676
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.22669914 -0.03305653
## sample estimates:
## mean in group 1 mean in group 2 
##       0.9199346       1.0498125
#効果量dの算出
library(effsize)
## 
## Attaching package: 'effsize'
## The following object is masked from 'package:psych':
## 
##     cohen.d
cohen.d(data$hsp_mean, data$gender)#HSP
## 
## Cohen's d
## 
## d estimate: -0.3555572 (small)
## 95 percent confidence interval:
##      lower      upper 
## -0.5447318 -0.1663826
cohen.d(data$eoe, data$gender)#EOE
## 
## Cohen's d
## 
## d estimate: -0.3659958 (small)
## 95 percent confidence interval:
##      lower      upper 
## -0.5552576 -0.1767340
cohen.d(data$lst, data$gender)#LST
## 
## Cohen's d
## 
## d estimate: -0.2311998 (small)
## 95 percent confidence interval:
##       lower       upper 
## -0.41952788 -0.04287168
cohen.d(data$aes, data$gender)#AES
## 
## Cohen's d
## 
## d estimate: -0.2106421 (small)
## 95 percent confidence interval:
##       lower       upper 
## -0.39886458 -0.02241967
cohen.d(data$brs_mean, data$gender)#BRS
## 
## Cohen's d
## 
## d estimate: 0.1488497 (negligible)
## 95 percent confidence interval:
##      lower      upper 
## -0.0391141  0.3368136
cohen.d(data$covid_mean, data$gender)#COVID
## 
## Cohen's d
## 
## d estimate: -0.2517916 (small)
## 95 percent confidence interval:
##       lower       upper 
## -0.44023529 -0.06334788

(2)階層的重回帰分析

2-1. 独立変数の中心化

data_c <- data %>% select("gender","age","hsp_mean","eoe","lst","aes","brs_mean","covid_mean") %>% drop_na()
names(data_c) #確認
## [1] "gender"     "age"        "hsp_mean"   "eoe"        "lst"       
## [6] "aes"        "brs_mean"   "covid_mean"
data_c$age <- data_c$age - mean(data_c$age) #ageの中心化
data_c$gender <- as.factor(data_c$gender)
data_c$hsp_mean <- data_c$hsp_mean - mean(data_c$hsp_mean) #hsp_meanの中心化
data_c$eoe <- data_c$eoe - mean(data_c$eoe) #eoeの中心化
data_c$lst <- data_c$lst - mean(data_c$lst) #lstの中心化
data_c$aes <- data_c$aes - mean(data_c$aes) #aesの中心化
data_c$brs_mean <- data_c$brs_mean - mean(data_c$brs_mean) #brs_meanの中心化
data_c$covid_mean <- data_c$covid_mean - mean(data_c$covid_mean) #covid_meanの中心化
names(data_c)
## [1] "gender"     "age"        "hsp_mean"   "eoe"        "lst"       
## [6] "aes"        "brs_mean"   "covid_mean"

2-2. HSP

##ステップ1:性別、年齢
model_s1 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender, data = data_c)
summary(model_s1)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender, 
##     data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04560 -0.36450 -0.07056  0.31551  2.13550 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    -0.06747    0.03624  -1.862   0.0633 .
## data_c$age     -0.02496    0.03008  -0.830   0.4071  
## data_c$gender2  0.12554    0.04956   2.533   0.0116 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.516 on 438 degrees of freedom
## Multiple R-squared:  0.01713,    Adjusted R-squared:  0.01264 
## F-statistic: 3.817 on 2 and 438 DF,  p-value: 0.02273
AIC(model_s1)
## [1] 672.916
BIC(model_s1)
## [1] 689.2722
##ステップ2:レジリエンス、HSP
model_s2 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender + data_c$brs_mean + data_c$hsp_mean, data = data_c)
summary(model_s2)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender + 
##     data_c$brs_mean + data_c$hsp_mean, data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13814 -0.34929 -0.06257  0.30947  2.00382 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     -0.04644    0.03533  -1.315  0.18936   
## data_c$age      -0.01273    0.02918  -0.436  0.66289   
## data_c$gender2   0.08641    0.04864   1.777  0.07631 . 
## data_c$brs_mean -0.08433    0.03177  -2.654  0.00823 **
## data_c$hsp_mean  0.07855    0.02467   3.184  0.00156 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4991 on 436 degrees of freedom
## Multiple R-squared:  0.08466,    Adjusted R-squared:  0.07626 
## F-statistic: 10.08 on 4 and 436 DF,  p-value: 8.209e-08
AIC(model_s2)
## [1] 645.5268
BIC(model_s2)
## [1] 670.061
anova(model_s1, model_s2) #R^2の増加量の検定
##ステップ3:レジリエンス*HSP交互作用
model_s3 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender + data_c$brs_mean + data_c$hsp_mean + data_c$brs_mean:data_c$hsp_mean, data = data_c)
summary(model_s3)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender + 
##     data_c$brs_mean + data_c$hsp_mean + data_c$brs_mean:data_c$hsp_mean, 
##     data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13919 -0.35032 -0.05829  0.30485  2.00321 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                     -0.04026    0.03690  -1.091  0.27588   
## data_c$age                      -0.01308    0.02921  -0.448  0.65452   
## data_c$gender2                   0.08684    0.04868   1.784  0.07511 . 
## data_c$brs_mean                 -0.08574    0.03188  -2.689  0.00744 **
## data_c$hsp_mean                  0.07818    0.02469   3.166  0.00166 **
## data_c$brs_mean:data_c$hsp_mean  0.01430    0.02444   0.585  0.55881   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4995 on 435 degrees of freedom
## Multiple R-squared:  0.08538,    Adjusted R-squared:  0.07486 
## F-statistic: 8.121 on 5 and 435 DF,  p-value: 2.426e-07
AIC(model_s3)
## [1] 647.1799
BIC(model_s3)
## [1] 675.8032
anova(model_s2, model_s3) #R^2の増加量の検定

2-3. EOE

##ステップ1:性別、年齢
model_s1 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender, data = data_c)
summary(model_s1)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender, 
##     data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04560 -0.36450 -0.07056  0.31551  2.13550 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    -0.06747    0.03624  -1.862   0.0633 .
## data_c$age     -0.02496    0.03008  -0.830   0.4071  
## data_c$gender2  0.12554    0.04956   2.533   0.0116 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.516 on 438 degrees of freedom
## Multiple R-squared:  0.01713,    Adjusted R-squared:  0.01264 
## F-statistic: 3.817 on 2 and 438 DF,  p-value: 0.02273
AIC(model_s1)
## [1] 672.916
BIC(model_s1)
## [1] 689.2722
##ステップ2:レジリエンス、EOE
model_s2 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender + data_c$brs_mean + data_c$eoe, data = data_c)
summary(model_s2)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender + 
##     data_c$brs_mean + data_c$eoe, data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13726 -0.35346 -0.05397  0.31346  2.04673 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     -0.04841    0.03548  -1.364   0.1732  
## data_c$age      -0.01116    0.02933  -0.381   0.7036  
## data_c$gender2   0.09007    0.04887   1.843   0.0660 .
## data_c$brs_mean -0.08214    0.03409  -2.410   0.0164 *
## data_c$eoe       0.06394    0.02474   2.585   0.0101 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.501 on 436 degrees of freedom
## Multiple R-squared:  0.07751,    Adjusted R-squared:  0.06904 
## F-statistic: 9.158 on 4 and 436 DF,  p-value: 4.119e-07
AIC(model_s2)
## [1] 648.9585
BIC(model_s2)
## [1] 673.4928
anova(model_s1, model_s2) #R^2の増加量の検定
##ステップ3:レジリエンス*EOE交互作用
model_s3 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender + data_c$brs_mean + data_c$hsp_mean + data_c$brs_mean:data_c$eoe, data = data_c)
summary(model_s3)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender + 
##     data_c$brs_mean + data_c$hsp_mean + data_c$brs_mean:data_c$eoe, 
##     data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13974 -0.35180 -0.06834  0.30091  1.99301 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                -0.02862    0.03777  -0.758  0.44890   
## data_c$age                 -0.01478    0.02920  -0.506  0.61303   
## data_c$gender2              0.08509    0.04860   1.751  0.08070 . 
## data_c$brs_mean            -0.08637    0.03178  -2.718  0.00683 **
## data_c$hsp_mean             0.07733    0.02466   3.136  0.00183 **
## data_c$brs_mean:data_c$eoe  0.02972    0.02241   1.326  0.18541   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4987 on 435 degrees of freedom
## Multiple R-squared:  0.08834,    Adjusted R-squared:  0.07787 
## F-statistic: 8.431 on 5 and 435 DF,  p-value: 1.256e-07
AIC(model_s3)
## [1] 645.7468
BIC(model_s3)
## [1] 674.3701
anova(model_s2, model_s3) #R^2の増加量の検定

2-4. LST

##ステップ1:性別、年齢
model_s1 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender, data = data_c)
summary(model_s1)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender, 
##     data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04560 -0.36450 -0.07056  0.31551  2.13550 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    -0.06747    0.03624  -1.862   0.0633 .
## data_c$age     -0.02496    0.03008  -0.830   0.4071  
## data_c$gender2  0.12554    0.04956   2.533   0.0116 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.516 on 438 degrees of freedom
## Multiple R-squared:  0.01713,    Adjusted R-squared:  0.01264 
## F-statistic: 3.817 on 2 and 438 DF,  p-value: 0.02273
AIC(model_s1)
## [1] 672.916
BIC(model_s1)
## [1] 689.2722
##ステップ2:レジリエンス、LST
model_s2 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender + data_c$brs_mean + data_c$lst, data = data_c)
summary(model_s2)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender + 
##     data_c$brs_mean + data_c$lst, data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13696 -0.34925 -0.06629  0.29666  1.96539 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.05181    0.03519  -1.472 0.141750    
## data_c$age      -0.01415    0.02919  -0.485 0.628120    
## data_c$gender2   0.09640    0.04828   1.997 0.046487 *  
## data_c$brs_mean -0.10180    0.02968  -3.430 0.000662 ***
## data_c$lst       0.04961    0.01595   3.111 0.001987 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4994 on 436 degrees of freedom
## Multiple R-squared:  0.08371,    Adjusted R-squared:  0.0753 
## F-statistic: 9.958 on 4 and 436 DF,  p-value: 1.018e-07
AIC(model_s2)
## [1] 645.9832
BIC(model_s2)
## [1] 670.5174
anova(model_s1, model_s2) #R^2の増加量の検定
##ステップ3:レジリエンス*LST交互作用
model_s3 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender + data_c$brs_mean + data_c$hsp_mean + data_c$brs_mean:data_c$lst, data = data_c)
summary(model_s3)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender + 
##     data_c$brs_mean + data_c$hsp_mean + data_c$brs_mean:data_c$lst, 
##     data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13831 -0.36439 -0.07142  0.31835  1.99766 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                -0.05216    0.03609  -1.446  0.14902   
## data_c$age                 -0.01313    0.02920  -0.450  0.65307   
## data_c$gender2              0.08620    0.04866   1.771  0.07718 . 
## data_c$brs_mean            -0.08273    0.03185  -2.598  0.00971 **
## data_c$hsp_mean             0.07854    0.02468   3.183  0.00156 **
## data_c$brs_mean:data_c$lst -0.01306    0.01662  -0.786  0.43241   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4993 on 435 degrees of freedom
## Multiple R-squared:  0.08595,    Adjusted R-squared:  0.07545 
## F-statistic: 8.181 on 5 and 435 DF,  p-value: 2.135e-07
AIC(model_s3)
## [1] 646.9012
BIC(model_s3)
## [1] 675.5245
anova(model_s2, model_s3) #R^2の増加量の検定

2-5. AES

##ステップ1:性別、年齢
model_s1 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender, data = data_c)
summary(model_s1)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender, 
##     data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04560 -0.36450 -0.07056  0.31551  2.13550 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    -0.06747    0.03624  -1.862   0.0633 .
## data_c$age     -0.02496    0.03008  -0.830   0.4071  
## data_c$gender2  0.12554    0.04956   2.533   0.0116 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.516 on 438 degrees of freedom
## Multiple R-squared:  0.01713,    Adjusted R-squared:  0.01264 
## F-statistic: 3.817 on 2 and 438 DF,  p-value: 0.02273
AIC(model_s1)
## [1] 672.916
BIC(model_s1)
## [1] 689.2722
##ステップ2:レジリエンス、AES
model_s2 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender + data_c$brs_mean + data_c$aes, data = data_c)
summary(model_s2)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender + 
##     data_c$brs_mean + data_c$aes, data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.10420 -0.34283 -0.09961  0.33365  2.06507 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.05585    0.03552  -1.573   0.1166    
## data_c$age      -0.01587    0.02945  -0.539   0.5901    
## data_c$gender2   0.10392    0.04873   2.133   0.0335 *  
## data_c$brs_mean -0.12798    0.02846  -4.497 8.86e-06 ***
## data_c$aes       0.02156    0.01569   1.374   0.1702    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5038 on 436 degrees of freedom
## Multiple R-squared:  0.06741,    Adjusted R-squared:  0.05885 
## F-statistic: 7.879 on 4 and 436 DF,  p-value: 3.878e-06
AIC(model_s2)
## [1] 653.76
BIC(model_s2)
## [1] 678.2943
anova(model_s1, model_s2) #R^2の増加量の検定
##ステップ3:レジリエンス*AES交互作用
model_s3 <- lm(data_c$covid_mean ~ data_c$age + data_c$gender + data_c$brs_mean + data_c$hsp_mean + data_c$brs_mean:data_c$aes, data = data_c)
summary(model_s3)
## 
## Call:
## lm(formula = data_c$covid_mean ~ data_c$age + data_c$gender + 
##     data_c$brs_mean + data_c$hsp_mean + data_c$brs_mean:data_c$aes, 
##     data = data_c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.14165 -0.34725 -0.05684  0.31576  2.00364 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                -0.04636    0.03534  -1.312  0.19032   
## data_c$age                 -0.01268    0.02919  -0.434  0.66414   
## data_c$gender2              0.08974    0.04885   1.837  0.06687 . 
## data_c$brs_mean            -0.08626    0.03188  -2.706  0.00708 **
## data_c$hsp_mean             0.07816    0.02468   3.166  0.00165 **
## data_c$brs_mean:data_c$aes  0.01396    0.01800   0.776  0.43838   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4993 on 435 degrees of freedom
## Multiple R-squared:  0.08592,    Adjusted R-squared:  0.07541 
## F-statistic: 8.178 on 5 and 435 DF,  p-value: 2.151e-07
AIC(model_s3)
## [1] 646.9173
BIC(model_s3)
## [1] 675.5406
anova(model_s2, model_s3) #R^2の増加量の検定

(3)媒介分析

3-1. HSP

library(lavaan)

#直接効果(c')のみのモデル
model <- "
      #コントロール
      hsp_mean ~ age + gender
      
      #直接効果
      covid_mean  ~ c*hsp_mean
"
fit <- sem(model, data = data, estimator = "ML")
summary(fit, standardized = TRUE, fit.measure = TRUE, ci = TRUE)
## lavaan 0.6-5 ended normally after 19 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                          5
##                                                       
##   Number of observations                           441
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                 3.532
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.171
## 
## Model Test Baseline Model:
## 
##   Test statistic                                47.073
##   Degrees of freedom                                 5
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.964
##   Tucker-Lewis Index (TLI)                       0.909
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)               -986.037
##   Loglikelihood unrestricted model (H1)             NA
##                                                       
##   Akaike (AIC)                                1982.075
##   Bayesian (BIC)                              2002.520
##   Sample-size adjusted Bayesian (BIC)         1986.652
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.042
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.112
##   P-value RMSEA <= 0.05                          0.471
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.028
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##   hsp_mean ~                                                            
##     age              -0.076    0.064   -1.195    0.232   -0.201    0.049
##     gender            0.376    0.105    3.591    0.000    0.171    0.581
##   covid_mean ~                                                          
##     hsp_mean   (c)    0.117    0.022    5.417    0.000    0.075    0.159
##    Std.lv  Std.all
##                   
##    -0.076   -0.056
##     0.376    0.169
##                   
##     0.117    0.250
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##    .hsp_mean          1.189    0.080   14.849    0.000    1.032    1.346
##    .covid_mean        0.252    0.017   14.849    0.000    0.219    0.286
##    Std.lv  Std.all
##     1.189    0.966
##     0.252    0.938
fitMeasures(fit)
##                npar                fmin               chisq                  df 
##               5.000               0.004               3.532               2.000 
##              pvalue      baseline.chisq         baseline.df     baseline.pvalue 
##               0.171              47.073               5.000               0.000 
##                 cfi                 tli                nnfi                 rfi 
##               0.964               0.909               0.909               0.812 
##                 nfi                pnfi                 ifi                 rni 
##               0.925               0.370               0.966               0.964 
##                logl   unrestricted.logl                 aic                 bic 
##            -986.037                  NA            1982.075            2002.520 
##              ntotal                bic2               rmsea      rmsea.ci.lower 
##             441.000            1986.652               0.042               0.000 
##      rmsea.ci.upper        rmsea.pvalue                 rmr          rmr_nomean 
##               0.112               0.471               0.008               0.008 
##                srmr        srmr_bentler srmr_bentler_nomean                crmr 
##               0.028               0.028               0.028               0.036 
##         crmr_nomean          srmr_mplus   srmr_mplus_nomean               cn_05 
##               0.036               0.028               0.028             749.127 
##               cn_01                 gfi                agfi                pgfi 
##            1151.053               0.992               0.960               0.198 
##                 mfi                ecvi 
##               0.998               0.031
#媒介変数を含めたモデル
set.seed(111)
model <- "
      #コントロール
      hsp_mean ~ age + gender
      
      #直接効果
      covid_mean  ~ c*hsp_mean
       
      #媒介
      brs_mean ~ a1*hsp_mean
      covid_mean  ~ b1*brs_mean
 
      #間接効果 (a*b)
      ab := a1*b1
       
      #全体の効果
      total := c + ab
"
fit <- sem(model, data = data, estimator = "ML", se = "bootstrap", bootstrap = 5000)
summary(fit, standardized = TRUE, fit.measure = TRUE, ci = TRUE)
## lavaan 0.6-5 ended normally after 21 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                          8
##                                                       
##   Number of observations                           441
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                 4.736
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.316
## 
## Model Test Baseline Model:
## 
##   Test statistic                               168.106
##   Degrees of freedom                                 9
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.995
##   Tucker-Lewis Index (TLI)                       0.990
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1480.879
##   Loglikelihood unrestricted model (H1)             NA
##                                                       
##   Akaike (AIC)                                2977.759
##   Bayesian (BIC)                              3010.471
##   Sample-size adjusted Bayesian (BIC)         2985.083
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.020
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.077
##   P-value RMSEA <= 0.05                          0.742
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.025
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             5000
##   Number of successful bootstrap draws            5000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##   hsp_mean ~                                                            
##     age              -0.076    0.060   -1.268    0.205   -0.190    0.043
##     gender            0.376    0.103    3.650    0.000    0.170    0.572
##   covid_mean ~                                                          
##     hsp_mean   (c)    0.086    0.025    3.394    0.001    0.037    0.137
##   brs_mean ~                                                            
##     hsp_mean  (a1)   -0.365    0.033  -11.042    0.000   -0.427   -0.299
##   covid_mean ~                                                          
##     brs_mean  (b1)   -0.084    0.033   -2.591    0.010   -0.149   -0.022
##    Std.lv  Std.all
##                   
##    -0.076   -0.056
##     0.376    0.169
##                   
##     0.086    0.184
##                   
##    -0.365   -0.475
##                   
##    -0.084   -0.139
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##    .hsp_mean          1.189    0.073   16.354    0.000    1.045    1.329
##    .covid_mean        0.248    0.019   12.918    0.000    0.211    0.286
##    .brs_mean          0.561    0.033   16.967    0.000    0.496    0.626
##    Std.lv  Std.all
##     1.189    0.966
##     0.248    0.923
##     0.561    0.774
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##     ab                0.031    0.012    2.488    0.013    0.008    0.056
##     total             0.117    0.022    5.424    0.000    0.076    0.159
##    Std.lv  Std.all
##     0.031    0.066
##     0.117    0.250
fitMeasures(fit)
##                npar                fmin               chisq                  df 
##               8.000               0.005               4.736               4.000 
##              pvalue      baseline.chisq         baseline.df     baseline.pvalue 
##               0.316             168.106               9.000               0.000 
##                 cfi                 tli                nnfi                 rfi 
##               0.995               0.990               0.990               0.937 
##                 nfi                pnfi                 ifi                 rni 
##               0.972               0.432               0.996               0.995 
##                logl   unrestricted.logl                 aic                 bic 
##           -1480.879                  NA            2977.759            3010.471 
##              ntotal                bic2               rmsea      rmsea.ci.lower 
##             441.000            2985.083               0.020               0.000 
##      rmsea.ci.upper        rmsea.pvalue                 rmr          rmr_nomean 
##               0.077               0.742               0.010               0.010 
##                srmr        srmr_bentler srmr_bentler_nomean                crmr 
##               0.025               0.025               0.025               0.031 
##         crmr_nomean          srmr_mplus   srmr_mplus_nomean               cn_05 
##               0.031               0.025               0.025             884.543 
##               cn_01                 gfi                agfi                pgfi 
##            1237.391               0.993               0.973               0.265 
##                 mfi                ecvi 
##               0.999               0.047

3-2. EOE

set.seed(222)
model <- "
      #コントロール
      eoe ~ age + gender
      
      #直接効果
      covid_mean  ~ c*eoe
       
      #媒介
      brs_mean ~ a1*eoe
      covid_mean  ~ b1*brs_mean
 
      #間接効果 (a*b)
      ab := a1*b1
       
      #全体の効果
      total := c + ab
"
fit <- sem(model, data = data, estimator = "ML", se = "bootstrap", bootstrap = 5000)
summary(fit, standardized = TRUE, fit.measure = TRUE, ci = TRUE)
## lavaan 0.6-5 ended normally after 22 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                          8
##                                                       
##   Number of observations                           441
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                 4.627
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.328
## 
## Model Test Baseline Model:
## 
##   Test statistic                               225.829
##   Degrees of freedom                                 9
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.997
##   Tucker-Lewis Index (TLI)                       0.993
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1483.026
##   Loglikelihood unrestricted model (H1)             NA
##                                                       
##   Akaike (AIC)                                2982.051
##   Bayesian (BIC)                              3014.764
##   Sample-size adjusted Bayesian (BIC)         2989.375
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.019
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.077
##   P-value RMSEA <= 0.05                          0.752
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.024
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             5000
##   Number of successful bootstrap draws            5000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##   eoe ~                                                                 
##     age              -0.120    0.063   -1.913    0.056   -0.250   -0.001
##     gender            0.409    0.112    3.644    0.000    0.193    0.625
##   covid_mean ~                                                          
##     eoe        (c)    0.072    0.026    2.777    0.005    0.021    0.122
##   brs_mean ~                                                            
##     eoe       (a1)   -0.406    0.028  -14.527    0.000   -0.461   -0.351
##   covid_mean ~                                                          
##     brs_mean  (b1)   -0.080    0.037   -2.172    0.030   -0.154   -0.009
##    Std.lv  Std.all
##                   
##    -0.120   -0.083
##     0.409    0.171
##                   
##     0.072    0.166
##                   
##    -0.406   -0.568
##                   
##    -0.080   -0.132
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##    .eoe               1.361    0.085   16.031    0.000    1.190    1.527
##    .covid_mean        0.250    0.020   12.502    0.000    0.212    0.291
##    .brs_mean          0.491    0.029   17.209    0.000    0.435    0.547
##    Std.lv  Std.all
##     1.361    0.961
##     0.250    0.930
##     0.491    0.677
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##     ab                0.033    0.015    2.145    0.032    0.003    0.064
##     total             0.105    0.019    5.436    0.000    0.067    0.142
##    Std.lv  Std.all
##     0.033    0.075
##     0.105    0.241
fitMeasures(fit)
##                npar                fmin               chisq                  df 
##               8.000               0.005               4.627               4.000 
##              pvalue      baseline.chisq         baseline.df     baseline.pvalue 
##               0.328             225.829               9.000               0.000 
##                 cfi                 tli                nnfi                 rfi 
##               0.997               0.993               0.993               0.954 
##                 nfi                pnfi                 ifi                 rni 
##               0.980               0.435               0.997               0.997 
##                logl   unrestricted.logl                 aic                 bic 
##           -1483.026                  NA            2982.051            3014.764 
##              ntotal                bic2               rmsea      rmsea.ci.lower 
##             441.000            2989.375               0.019               0.000 
##      rmsea.ci.upper        rmsea.pvalue                 rmr          rmr_nomean 
##               0.077               0.752               0.008               0.008 
##                srmr        srmr_bentler srmr_bentler_nomean                crmr 
##               0.024               0.024               0.024               0.029 
##         crmr_nomean          srmr_mplus   srmr_mplus_nomean               cn_05 
##               0.029               0.024               0.024             905.242 
##               cn_01                 gfi                agfi                pgfi 
##            1266.355               0.993               0.974               0.265 
##                 mfi                ecvi 
##               0.999               0.047

3-3. LST

set.seed(333)
model <- "
      #コントロール
      lst ~ age + gender
      
      #直接効果
      covid_mean  ~ c*lst
       
      #媒介
      brs_mean ~ a1*lst
      covid_mean  ~ b1*brs_mean
 
      #間接効果 (a*b)
      ab := a1*b1
       
      #全体の効果
      total := c + ab
"
fit <- sem(model, data = data, estimator = "ML", se = "bootstrap", bootstrap = 5000)
summary(fit, standardized = TRUE, fit.measure = TRUE, ci = TRUE)
## lavaan 0.6-5 ended normally after 25 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                          8
##                                                       
##   Number of observations                           441
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                 6.947
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.139
## 
## Model Test Baseline Model:
## 
##   Test statistic                                98.501
##   Degrees of freedom                                 9
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.967
##   Tucker-Lewis Index (TLI)                       0.926
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1674.863
##   Loglikelihood unrestricted model (H1)             NA
##                                                       
##   Akaike (AIC)                                3365.725
##   Bayesian (BIC)                              3398.438
##   Sample-size adjusted Bayesian (BIC)         3373.049
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.041
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.090
##   P-value RMSEA <= 0.05                          0.547
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.034
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             5000
##   Number of successful bootstrap draws            5000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##   lst ~                                                                 
##     age              -0.066    0.088   -0.746    0.456   -0.226    0.119
##     gender            0.354    0.153    2.317    0.021    0.052    0.652
##   covid_mean ~                                                          
##     lst        (c)    0.053    0.017    3.077    0.002    0.019    0.086
##   brs_mean ~                                                            
##     lst       (a1)   -0.177    0.026   -6.889    0.000   -0.228   -0.127
##   covid_mean ~                                                          
##     brs_mean  (b1)   -0.105    0.028   -3.697    0.000   -0.160   -0.047
##    Std.lv  Std.all
##                   
##    -0.066   -0.034
##     0.354    0.111
##                   
##     0.053    0.162
##                   
##    -0.177   -0.331
##                   
##    -0.105   -0.172
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##    .lst               2.485    0.122   20.294    0.000    2.224    2.702
##    .covid_mean        0.249    0.019   13.343    0.000    0.212    0.285
##    .brs_mean          0.645    0.037   17.482    0.000    0.569    0.717
##    Std.lv  Std.all
##     2.485    0.986
##     0.249    0.926
##     0.645    0.891
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##     ab                0.019    0.006    3.160    0.002    0.008    0.031
##     total             0.071    0.016    4.394    0.000    0.040    0.103
##    Std.lv  Std.all
##     0.019    0.057
##     0.071    0.219
fitMeasures(fit)
##                npar                fmin               chisq                  df 
##               8.000               0.008               6.947               4.000 
##              pvalue      baseline.chisq         baseline.df     baseline.pvalue 
##               0.139              98.501               9.000               0.000 
##                 cfi                 tli                nnfi                 rfi 
##               0.967               0.926               0.926               0.841 
##                 nfi                pnfi                 ifi                 rni 
##               0.929               0.413               0.969               0.967 
##                logl   unrestricted.logl                 aic                 bic 
##           -1674.863                  NA            3365.725            3398.438 
##              ntotal                bic2               rmsea      rmsea.ci.lower 
##             441.000            3373.049               0.041               0.000 
##      rmsea.ci.upper        rmsea.pvalue                 rmr          rmr_nomean 
##               0.090               0.547               0.015               0.015 
##                srmr        srmr_bentler srmr_bentler_nomean                crmr 
##               0.034               0.034               0.034               0.041 
##         crmr_nomean          srmr_mplus   srmr_mplus_nomean               cn_05 
##               0.041               0.034               0.034             603.290 
##               cn_01                 gfi                agfi                pgfi 
##             843.817               0.990               0.961               0.264 
##                 mfi                ecvi 
##               0.997               0.052

3-4. AES

set.seed(444)
model <- "
      #コントロール
      aes ~ age + gender
      
      #直接効果
      covid_mean  ~ c*aes
       
      #媒介
      brs_mean ~ a1*aes
      covid_mean  ~ b1*brs_mean
 
      #間接効果 (a*b)
      ab := a1*b1
       
      #全体の効果
      total := c + ab
"
fit <- sem(model, data = data, estimator = "ML", se = "bootstrap", bootstrap = 5000)
summary(fit, standardized = TRUE, fit.measure = TRUE, ci = TRUE)
## lavaan 0.6-5 ended normally after 23 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                          8
##                                                       
##   Number of observations                           441
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                 9.315
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.054
## 
## Model Test Baseline Model:
## 
##   Test statistic                                44.427
##   Degrees of freedom                                 9
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.850
##   Tucker-Lewis Index (TLI)                       0.662
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1690.866
##   Loglikelihood unrestricted model (H1)             NA
##                                                       
##   Akaike (AIC)                                3397.732
##   Bayesian (BIC)                              3430.444
##   Sample-size adjusted Bayesian (BIC)         3405.056
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.055
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.102
##   P-value RMSEA <= 0.05                          0.363
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.042
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             5000
##   Number of successful bootstrap draws            5000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##   aes ~                                                                 
##     age               0.019    0.091    0.212    0.832   -0.158    0.193
##     gender            0.328    0.150    2.184    0.029    0.033    0.621
##   covid_mean ~                                                          
##     aes        (c)    0.025    0.017    1.437    0.151   -0.009    0.059
##   brs_mean ~                                                            
##     aes       (a1)   -0.056    0.026   -2.123    0.034   -0.108   -0.004
##   covid_mean ~                                                          
##     brs_mean  (b1)   -0.133    0.028   -4.799    0.000   -0.186   -0.079
##    Std.lv  Std.all
##                   
##     0.019    0.010
##     0.328    0.106
##                   
##     0.025    0.074
##                   
##    -0.056   -0.102
##                   
##    -0.133   -0.218
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##    .aes               2.358    0.125   18.922    0.000    2.097    2.592
##    .covid_mean        0.254    0.019   13.194    0.000    0.217    0.292
##    .brs_mean          0.717    0.039   18.166    0.000    0.639    0.793
##    Std.lv  Std.all
##     2.358    0.989
##     0.254    0.944
##     0.717    0.990
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
##     ab                0.007    0.004    1.876    0.061    0.000    0.016
##     total             0.032    0.017    1.872    0.061   -0.001    0.067
##    Std.lv  Std.all
##     0.007    0.022
##     0.032    0.096
fitMeasures(fit)
##                npar                fmin               chisq                  df 
##               8.000               0.011               9.315               4.000 
##              pvalue      baseline.chisq         baseline.df     baseline.pvalue 
##               0.054              44.427               9.000               0.000 
##                 cfi                 tli                nnfi                 rfi 
##               0.850               0.662               0.662               0.528 
##                 nfi                pnfi                 ifi                 rni 
##               0.790               0.351               0.869               0.850 
##                logl   unrestricted.logl                 aic                 bic 
##           -1690.866                  NA            3397.732            3430.444 
##              ntotal                bic2               rmsea      rmsea.ci.lower 
##             441.000            3405.056               0.055               0.000 
##      rmsea.ci.upper        rmsea.pvalue                 rmr          rmr_nomean 
##               0.102               0.363               0.018               0.018 
##                srmr        srmr_bentler srmr_bentler_nomean                crmr 
##               0.042               0.042               0.042               0.051 
##         crmr_nomean          srmr_mplus   srmr_mplus_nomean               cn_05 
##               0.051               0.042               0.042             450.165 
##               cn_01                 gfi                agfi                pgfi 
##             629.541               0.986               0.948               0.263 
##                 mfi                ecvi 
##               0.994               0.057